Soil phosphorus cycling microbial functional genes of monoculture and mixed plantations of native tree species in subtropical China

Background Transforming coniferous plantation into broadleaved or mixed broadleaved-coniferous plantations is the tendency of forest management strategies in subtropical China. However, the effects of this conversion on soil phosphorus (P) cycling microbial functional genes are still unknown. Methods Soil samples were collected from 0–20, 20–40, and 40–60 cm (topsoil, middle layer, and subsoil, respectively) under coniferous Pinus massoniana (PM), broadleaved Erythrophleum fordii (EF), and their mixed (PM/EF) plantation in subtropical China. Used metagenomic sequencing to examine the alterations of relative abundances and molecular ecological network structure of soil P-cycling functional genes after the conversion of plantations. Results The composition of P-cycling genes in the topsoil of PM stand was significantly different from that of PM/EF and EF stands (p < 0.05), and total phosphorus (TP) was the main factor causing this difference. After transforming PM plantation into EF plantation, the relative abundances of P solubilization and mineralization genes significantly increased in the topsoil and middle layer with the decrease of soil TP content. The abundances of P-starvation response regulation genes also significantly increased in the subsoil (p < 0.05), which may have been influenced by soil organic carbon (SOC). The dominant genes in all soil layers under three plantations were phoR, glpP, gcd, ppk, and ppx. Transforming PM into EF plantation apparently increased gcd abundance in the topsoil (p < 0.05), with TP and NO3−-N being the main influencing factors. After transforming PM into PM/EF plantations, the molecular ecological network structure of P-cycling genes was more complex; moreover, the key genes in the network were modified with the transformation of PM plantation. Conclusion Transforming PM into EF plantation mainly improved the phosphate solubilizing potential of microorganisms at topsoil, while transforming PM into PM/EF plantation may have enhanced structural stability of microbial P-cycling genes react to environmental changes.


Introduction
Phosphorus (P) is a key macronutrient found in all living organisms and plays an essential role in plant growth and development (Reinhard et al., 2017).Although P is relatively abundant in soil, 95-99% of it is immobilized and thus can hardly be absorbed or utilized by plants (Lucero et al., 2021), making it a limiting nutrient for both soil microorganisms and plants (Wang et al., 2021).Phosphorus limitation in terrestrial ecosystems is considered a major issue that needs to be urgently addressed for ecosystem management and restoration (Liang et al., 2020).
Microorganisms exert important effects on soil P cycling and regulate its availability (Zaidi et al., 2009).The microbial P-transformation processes can be mostly regulated by 3 groups of microbial genes, namely, genes involved in inorganic P-solubilization and organic P-mineralization, P-uptake and transport, and P-starvation response regulation (Bergkemper et al., 2016).The genes related to P-starvation response regulation (phoB, phoR, and phoU) allow microorganisms to utilize P (Dai et al., 2020) and they also control the expression of genes related to organic P-mineralization (phoD) and P-uptake and transport (pst) (Hsieh and Wanner, 2010).Microorganisms have efficient P uptake systems, and thus, they compete with other soil biota for available P (Bergkemper et al., 2016;Wu et al., 2022).Genes encoding low-affinity inorganic phosphate transporter (Pit) and high-affinity phosphate transporter (Pst) play a significant role in promoting microbial phosphorus uptake and transport (Bergkemper et al., 2016).Microorganisms that carry genes related to inorganic P-solubilization and organic P-mineralization have been found to produce enzymes involved in mineralizing organic P or organic anions and solubilizing inorganic P (Dai et al., 2020).The Gene gcd is an important molecular biomarker for soil P-solubilizing microbes (Rawat et al., 2020), and phoD is a common gene responsible for mineralizing organic P (Hu et al., 2020) and is usually analyzed within the soil metagenome (Tan et al., 2012).
In recent years, soil microbial phosphorus cycling functional genes have attracted increasing attention, with the primary focus on agricultural ecosystems (Liu J. et al., 2018;Dai et al., 2020;Siles et al., 2022).Soil physicochemical properties are affected by various factors such as litter and root exudates from different vegetation types, and changes in their amounts drive the structure and functional responses of the soil microbial community (Bardgett and van der Putten, 2014).At the same time, soil physicochemical properties vary significantly with soil depth (Cade-Menun et al., 2017), causing alterations in soil microbial community abundance and structure at different soil depths (Becerra et al., 2014).Nonetheless, it remains largely unclear whether the abundance and structure of soil phosphorus cycling microbial functional genes vary among different stand types with soil layer.In addition, soil phosphorus cycling occurs as the result of an interaction between different microbial functional genes (Kelliher et al., 2018), and clarifying this interaction is the basis for understanding microbial community stability (Zhou et al., 2010).Molecular ecological network analysis used to characterize microbial functional genes provides a robust method to decipher the potential interactions of functional genes in complex microbial communities (Zhou et al., 2010;Liang et al., 2016).It has also been reported that the molecular ecological network structure of soil phosphorus cycling functional genes in natural evergreen-deciduous broadleaved mixed forests is more complex than that in natural evergreen broadleaved and deciduous broadleaved forests, indicating that soil phosphorus cycling microbial functional genes and their associated microbial communities in natural broadleaved mixed forests are more stable (Cheng et al., 2021).The key nodes in molecular ecological networks have an important and exclusive role in maintaining community stability, while the removal of these nodes from the networks dramatically changes microbial structure and functions (Banerjee et al., 2018).Therefore, identifying key microbial functional genes can shed light on the responses of microbial communities to environmental changes.However, the molecular ecological network structure of phosphorus cycling microbial functional genes in different planted forests is still unclear.
The establishment of native broadleaved plantations and coniferous-broadleaved mixed plantations is gradually becoming the most promising forest management approach to replace coniferous plantations distributed across large areas of the world (Liu S. et al., 2018;Wang et al., 2018) and also to improve soil fertility and forest productivity (Gillespie et al., 2021).Multiple generations of coniferous plantations can reduce litter decomposition rate, nutrient turnover rate, soil fertility, and net primary productivity (Schall and Ammer, 2012), while native broadleaved tree species have higher litter decomposition rate and soil fertility (You et al., 2018a), and the development of mixed-species plantations can promote soil nutrient turnover and improve woodland productivity (Handa et al., 2014).Phosphorus has a low availability in subtropical forest soils (You et al., 2020), and it is an important nutrient limiting plant development (Bergkemper et al., 2016).Therefore, understanding the effects of different planted forest management models on soil microbial phosphorus cycling functional genes can provide a reference for the selection of tree species and a scientific basis for the tree species distribution in subtropical regions.In this paper, the typical native coniferous Pinus massoniana (PM) plantation, broadleaved Erythrophleum fordii (EF) plantation, and PM/EF mixed plantation were studied in the subtropical region of China.Based on the metagenomic sequencing data of various soil depths (0-20, 20-40, and 40-60 cm) under all plantation stands, this study aims to investigate: (1) the alterations of the composition structure and relative abundance of soil microbial phosphorus cycling functional genes and the dominant soil physicochemical factors affected by transforming the coniferous plantation into mixed broadleavedconiferous or broadleaved plantations; (2) the response profiles of the molecular ecological network structure of soil microbial phosphorus cycling functional genes to the conversion of the coniferous plantation to broadleaved or mixed broadleaved-coniferous plantations and also the potential keys genes in the network under three plantation types.a typical subtropical monsoon climate and two distinguishable dry and wet seasons (the rainy season occurring from April to September and the dry season from October to March of the next year).The annual mean rainfall is 1,350 mm and relative humidity is above 80%.The landform types are mainly low mountains and hills, and the soil type is lateritic soil (Luo et al., 2014), which is the same as oxisol in the USDA Soil Taxonomy.
Three typical stands of PM, EF, and PM/EF were selected as research objects in this study (Table 1).These three stands were replanted on clear-cut sites of the PM stand in 2006, with a planting density of 2,500 plants/hm2 and 3:1 as the ratio of PM/EF stands.In January 2021, three sample plots (20 m × 20 m) were set up for each plantation, and plots were spaced at least 20 m apart to avoid spatial autocorrelation.In each sample plot, soil samples at the depths 0-20, 20-40, and 40-60 cm (topsoil, middle layer, and subsoil, respectively) (He et al., 2013;You et al., 2018b) were obtained using the soil drill (with an inner diameter of 5 cm) according to the 5-point "S-type" sampling method.Soil samples taken from the same layer were blended into one mixed sample.A total of 27 soil samples (3 stands × 3 plots × 3 soil layers) were obtained, placed into the sterile sampling bag, preserved on the biological ice pack, and quickly transferred to the laboratory.Each freshly collected mixed soil sample was divided into three parts, with one part being preserved in the refrigerator at −20°C following the filtration with the 2-mm sieve and then used for metagenomic sequencing.The second part was sieved with the 2-mm mesh and preserved in the refrigerator at 4°C to determine nitrate nitrogen (NO 3 − -N) and ammonium nitrogen (NH 4 + -N).The third part was sieved with the 0.25-mm mesh after natural air drying and then used for the determination of other soil physicochemical properties.

Soil physicochemical parameters
We measured soil pH at the soil: water ratio of 1: 2.5 (w/v) with the pH meter (PHS-3C, Shanghai Jinhuan Instrument Co., Ltd., Shanghai, China).The potassium dichromate oxidation-external heating method was employed for determining soil organic carbon (SOC).The melt-molybdenum, antimony, and scandium colorimetry were used to measure total phosphorus (TP), whereas diacid (HCl-H 2 SO 4 ) was used to extract available phosphorus (AP) which was then determined by enzymolysis (INFINITE M200 PRO, TECAN, Switzerland).The above-mentioned methods for the determination of soil physicochemical properties are described in detail by Bao (2000).Auto Analyzer (AutoAnalyzer3, SEAL, Germany) was used to measure total nitrogen (TN) content in the soil after de-boiling with H 2 SO 4 .KCl was used to extract NH 4 + -N and NO 3 − -N, which were then measured with Auto Analyzer (Wu et al., 2019).TOC analyzer (Multi N/C 3100, Analytik Jena AG, Germany) was used for measuring the dissolved organic carbon (DOC) content.The carbon-to-nitrogen ratio (C/N) in the soil is defined as the ratio of the mass of soil organic carbon to the mass of total nitrogen.

DNA isolation and metagenomic sequencing
The FastDNA ® Spin Kit for Soil was used for DNA extraction.DNA concentration, purity, and integrity were evaluated by electrophoresis on a 1% agarose gel.The genomic DNA was segmented into 400-bp fragments using the Covaris M220 instrument.The NEXTFLEX Rapid DNA-Seq Kit was adopted for paired-end (PE) library construction.The metagenomic sequencing was performed on the Illumina NovaSeq sequencing platform at Shanghai Meiji Biotechnology Co., Ltd.Besides, we applied the fastp tool v0.20.0 to quality control checks on original sequence data and optimize sequences for metagenomic analysis.Megahit software v1.1.2was used to splice and assemble the optimized sequences, and overlapping groups of ≥300 bp from the splicing analysis were selected as the final assembly results for the next steps which were the gene prediction and annotation.MetaGene software was used to predict overlapping groups obtained from splicing analysis.We chose the genes with the nucleic acid length of ≥100 bp, translated them into amino-acid sequences, and finally obtained gene prediction statistics for every sample.Then, CD-HIT v4.6.1 was applied for clustering the gene sequences predicted in every sample (sequence similarity ≥90% and coverage ≥90%), while genes with the greatest length per class were chosen as the typical sequences for constructing the non-redundant gene set.SOAPaligner v2.21 was used to align high-quality sequences obtained from quality inspection of every sample and the non-redundant gene set (sequence similarity ≥95%).The abundance of each functional gene in the corresponding soil sample was calculated, and then the ratio of the abundance of each functional gene to the sum of the abundance of all functional genes in the sample was taken as the relative abundance of a single functional gene in the sample.The metagenomic sequencing data were imported into the National Microbiology Data Center1 (NMDC10018441).
Using BLASTP (BLAST Version 2.2.28+),2 sequences of redundant genes were compared with those of KEGG database, 3 and the e-value of BLAST comparison parameters is set to 1e-5.According to the comparison results, KOBAS 2.0 (Xie et al., 2011) (KEGG Orthology Based Annotation System) was used for functional annotation.A total of 55 microbial phosphorus cycling functional genes obtained from all soil samples were grouped into three main categories (phosphorusstarvation response regulation, phosphorus absorption and transport, and inorganic phosphorus solubilization and organic phosphate mineralization) (Table 2).

Statistical analysis
Effects of stand type, soil layer, and their interaction on the relative abundances of phosphorus cycling function genes and soil physicochemical properties were determined by the two-way analysis of variance (ANOVA), and Duncan's multiple range test (MRT) was used for multiple comparisons of different treatment groups.Before performing ANOVA, all data were subjected to Levene's test for the homogeneity of variance and Shapiro -Wilk normality test, and data that did not conform to the assumptions of ANOVA were transformed to log10.The calculations were carried out by using SPSS 26.0 (IBM SPSS Inc., Chicago, IL, USA).According to the Bray-Curtis distance calculated from the relative abundances of phosphorus cycling functional genes in each soil sample, principal coordinate analysis (PCoA), permutational multivariate analysis of variance (PERMANOVA), and redundancy analysis (RDA) were performed to investigate the differences in the composition structure of soil phosphorus cycling microbial functional genes in different stands and the main soil physicochemical factors leading to these differences.Collinear variables whose variance inflation factor (VIF) was >10 were excluded from the RDA model (Liu J. et al., 2018).The calculation was performed by the "cmdscale ()" function, "adonis ()" function, and "dbrda ()" function of the "vegan" package in R software v 4.0.5 (Dixon, 2003).
Pearson's correlation was employed to analyze the relationships between soil physicochemical properties and relative abundances of phosphorus cycling functional genes.The correlation heatmap was created by the "corr.Test ()" function of the psych package in R software v 4.0.5.Using the "lavaan" package in R software, we performed structural equation modeling (SEM) to quantify the effects of stand types on soil phosphorus cycling functional gene groups.Topological parameters of the molecular ecological network were calculated in the MENA platform4 for P-cycling microbial functional genes, and network visualization was carried out using Cytoscape 3.8.0software. 5The nodes of each network were classified into four categories based on their topological roles: (1) peripherals (Zi < 2.5 and Pi < 0.62); (2) connectors (Zi < 2.5 and Pi ≥ 0.62); (3) module hubs (Zi ≥ 2.5 and Pi < 0.62); and (4) network hubs (Zi ≥ 2.5 and Pi ≥ 0.62).Among them, the last three are regarded as the key nodes in molecular ecological networks (Olesen et al., 2007).

Soil properties
According to the two-way ANOVA, stand type only dramatically affected soil TP and TN, while soil layer remarkably affected SOC, DOC, TP, and C/N, and the interaction between stand type and soil layer apparently affected only soil TP (p < 0.05) (Table 3).The soil TP contents varied among PM, EF, and PM/EF stands with soil layer.In the topsoil (0-20 cm), the TP content of the PM stand remarkably increased compared to that of PM/EF and EF stands (p < 0.05).In the middle layer and subsoil (20-40 and 40-60 cm, respectively), the TP content of the PM stand had no significant difference from that of the PM/EF stand, but remarkably increased compared to that of the EF stand.In addition, TN and NO 3 − -N contents in the subsoil of the PM plantation were significantly lower than those of the PM/EF plantation.Differences in the other 6 soil physicochemical parameters, including pH, SOC, DOC, NO 3 − -N, AP, and C/N in all soil layers were not significant among the three plantation stands.

Composition structure of microbial genes involved in P transformation
In the topsoil, principal coordinate analysis (PCoA) showed that the first axis (PCoA1) could explain 43.97% of the variation of the structure of soil phosphorus cycling microbial functional genes, while the second axis (PCoA2) accounted for 18.34% of this variation, and the total interpretation rate of the first two axes reached 62.31% (Figure 1A).The soil samples from the PM plantation were located on the right side of PCoA1, whereas the EF plantation soil samples were located on the left side of PCoA1, and the soil samples of the PM/EF plantation were concentrated below PCoA2.PERMANOVA analysis further demonstrated that the structure of soil phosphorus cycling genes in the PM plantation was significantly different from that in both EF and PM/EF plantations (p < 0.05).RDA analysis showed that the first two axes explained 66.13% (RDA1:51.52%and RDA2:14.61%) of the impact of soil physicochemical factors on the composition structure of phosphorus cycling functional genes (Figure 1B).Monte Carlo test revealed that soil TP was the main physicochemical factor that caused the significant difference in the composition structure of soil phosphorus cycling microbial functional genes among PM, EF, and PM/EF plantations (p < 0.05).
In the middle layer, the cumulative interpretation rate of the first two axes of PCoA was 53.42%, and the interpretation rates of PCoA1 and PCoA2 were 36.28 and 17.14%, respectively (Figure 1C); however, 5 http://manual.cytoscape.org/en/stable/PERMANOVA analysis showed no significant difference among PM, EF, and PM/EF stands.According to the RDA analysis, the first two axes accounted for 41.62 and 29.94%, respectively, of the impacts of soil physical and chemical factors on the structure of soil phosphorus cycling microbial functional genes, explaining the total variation of 71.56% (Figure 1D).Monte Carlo test results showed that the differences in the structure of phosphorus cycling microbial genes among soil samples were not caused by the selected physical and chemical properties.
In the subsoil, however, the total interpretation rate of the first two axes of PCoA was 65.86% (PCoA1:49.26%and PCoA2:16.60%)(Figure 1E).Based on PERMANOVA results, PM, EF, and PM/EF plantation stands were not significantly different.Soil physical and chemical factors accumulated in the first two axes of RDA could explain 66.98% of variation in the composition structure of soil phosphorus cycling microbial functional genes (RDA1: 48.30% and RDA2: 18.68%) (Figure 1F).As revealed by Monte Carlo test results, soil pH were the major factors driving differences in the composition structure of phosphorus cycling functional genes in different soil samples (p < 0.05).

Categories of phosphorus cycling microbial functional genes
Based on the two-way ANOVA, stand types and soil layers dramatically affected the relative abundance of genes related to P-starvation response regulation, inorganic P-solubilization, and organic P-mineralization (p < 0.05), while relative abundances of genes involved in P-uptake and transportation were not affected by stand, soil layer, and their interaction (Table 4).The genes involved in P-starvation response regulation were not significantly different with regard to their relative abundances in the soil surface and middle layer, whereas, in the subsoil, the values remarkably decreased in the PM stand compared to those in the EF stand (p < 0.05) (Figure 2A).For genes involved in P-uptake and transportation, similarly, relative abundances did not show significant differences across three stands throughout the soil profile (Figure 2B).For genes involved in inorganic P-solubilization and organic P-mineralization, in contrast, relative abundances in the topsoil were significantly lower in the PM stand than in the EF stand.In the middle layer, the PM stand exhibited remarkably decreased values compared to EF and PM/EF stands, and no significant difference was observed among the three stands in the subsoil (Figure 2C).According to the structural equation modeling (SEM) results (Figure 3), p-value of Chi-square > 0.05, the comparative fit index (CFI) and goodness-of-fit index (GFI) were both >0.9, standardized root mean square residual (SRMR) < 0.09, indicating a good model fit (Hooper et al., 2008).There were no significant effects of stand types on genes involved in P-starvation response regulation and P-uptake and transportation, but they had direct significant positive effects on genes involved in inorganic P-solubilization and organic P-mineralization (p < 0.05).At the same time, stand types exerted indirect negative impacts on genes involved in inorganic P-solubilization and organic P-mineralization through the soil TP content.

Species of phosphorus cycling microbial functional genes
Predominant functional genes in every soil layer for 3 stands were phoR, glpP, gcd, ppk, and ppx, with a relative abundance of >0.5‰ (Figure 4).The two-way ANOVA showed that stand types and soil layers remarkably affected the relative abundances of phoR, gcd, and phnW (p < 0.05), while ugpQ and phoD were significantly affected only by soil layers, and phnX was significantly affected only by stand types (Table 4).In addition, at the topsoil, the relative abundances of gcd, ugpQ, phnW, and phnX in the PM plantation apparently decreased compared to those in the EF plantation (p < 0.05) (Figure 4A); at the middle layer, similarly, the relative abundance of phoD in the PM plantation remarkably decreased compared with that of EF and PM/ EF plantations (Figure 4B), and at the subsoil, the relative abundance of phoR in the PM plantation evidently decreased compared with that in the EF plantation (Figure 4C).
As revealed by Pearson's correlation analysis results, gcd showed an obvious negative relationship with TP but a positive relationship with NO 3 − -N, whereas phnW exhibited a negative relationship with AP, TP, and pH, and likewise, phnX was negatively correlated with TP at the topsoil.At the middle layer, phoD was negatively correlated with TP.However, there was no correlation between phoR and evaluated soil physicochemical properties.It can be seen that almost all genes with significant differences in their relative abundance among the three stands had significant negative correlations with TP (Figure 5).

Molecular ecological network analysis of microbial genes related to P transformation
The molecular ecological network analysis of phosphorus cycling microbial functional genes within soil profiles (0-60 cm) under all plantation stands revealed that 41, 34, and 41 gene nodes and 112, 103, and 179 lines were retained for PM, EF, and PM/EF stands, respectively, to construct molecular ecological networks of phosphorus cycling functional genes (Table 5).The connecting lines in three molecular ecological networks represent mainly positive interactions between nodes (Figure 6).
Topological analysis of molecular ecological networks of soil phosphorus cycling microbial functional genes in 3 studied stands showed that the complexity of the network (average degree), the degree of community clustering (mean clustering coefficient), the distance between nodes (mean path length), and level of community organization (modularity) increased after the transformation of the PM stand into the PM/EF stand.However, after the transformation of the PM plantation stand into the EF stand, only the mean clustering coefficient of the network increased, while mean path length, modularity, and average connectivity decreased.
In a molecular ecological network, different genes (nodes) in a microbial community play distinct roles (Fuhrman and Steele, 2008).To determine how plantation conversion affected the topological role of nodes in a network, within-and among-module connectivity (Zi and Pi, respectively) in three networks were visualized (Figure 7).There was only one module hub (ugpC) in the network for the PM plantation, while the network for the EF plantation was found to have four connectors (phoU, pstA, gcd, and     6).

Discussion
4.1 Effects of stand types on the composition structure of genes related to P transformation Differences in the biological characteristics of tree species in forest soils lead to differences in the quantity and quality of litter, thereby affecting soil microbial community structure within the forest floor (Richter et al., 2018).Changes in the composition of microbial functional genes are often associated with alterations in plant and soil characteristics in ecosystems (Hu et al., 2019).Soil TP content is an important factor that affects the composition structure of soil phosphorus cycling microbial functional genes (Yu et al., 2021;Wang et al., 2022).According to our results, there were significant differences among the PM, EF, and PM/EF plantations in terms of the composition structure of phosphorus-cycling microbial functional genes at the topsoil.The main reason for this difference is that the TP content in the PM stand (0.40 g‧kg −1 ) remarkably increased compared to that of the EF stand (0.28 g‧kg −1 ) and the PM/EF stand (0.33 g‧kg −1 ).It is worth noting that in this study, functional gene composition was not significantly different among three stands in the middle layer and subsoil, which indicates that the soil layer classification provides a framework that can deepen our understanding of how stand types affect the composition structure of soil phosphorus cycling microbial functional genes.

Effects of stand types on relative abundances of genes related to P transformation
According to our results, the relative abundance of P-starvation response regulation gene groups in the subsoil significantly decreased in the PM stand compared with that in the EF stand, which is probably attributed to the higher SOC content in this soil layer of the PM stand (16.07 g‧kg −1 ) than that in the EF stand (14.87 g‧kg −1 ).The SOC content showed a negative relationship with the relative abundance of the P-starvation response regulation genes (Li et al., 2022).In addition, concerning the genes involved in inorganic P-solubilization and organic P-mineralization, relative abundances in the topsoil and the middle layer apparently decreased in the PM stand compared to that in the EF stand.Structural equation modeling revealed that stand types could affect soil TP content, which in turn, had a significant negative effect on relative abundances of genes related to inorganic P-solubilization and organic P-mineralization.This is consistent with the results of previous studies indicating that under low-P conditions, genes associated with inorganic P-solubilization and organic P-mineralization enhanced phosphorus availability by encoding phosphatase or releasing organic acids, which can explain the negative correlation of soil TP content with the abundance of these genes (Luo et al., 2020).Additionally, the TP content markedly increased in the topsoil and the middle layer of the PM stand (0.4 and 0.3 g‧kg −1 , respectively) compared to that of the EF stand (0.28 and 0.23 g‧kg −1 , respectively), demonstrating it as the main factor leading to the above-mentioned difference.
Dominant functional genes associated with soil phosphorus cycling varied in different stand types (Cheng et al., 2022;Wu et al., FIGURE 3 Structural equation model (SEM) illustrating the effects of stand types on groups of phosphorus cycling microbial functional genes; Arrow width indicates the strength of the relationship.Blue and red colors represent positive and negative relationships, respectively.Solid and dashed arrows represent significant and non-significant relationships among different variables, respectively (*p < 0.05, **p < 0.01, and ***p < 0.001).The values adjacent to arrows indicate path coefficients.r 2 : percentage of the variance in each dependent variable.Qin et al. 10.3389/fmicb.2024.1419645Frontiers in Microbiology 11 frontiersin.org2022).Among the 55 functional genes evaluated in this study, the dominant ones in all soil layers of 3 plantations were phoR, glpP, gcd, ppk, and ppx, indicating that the dominant functional genes involved in phosphorus cycling in all soil layers (0-60 cm) were not modified after the transformation of the PM plantation into EF and PM/EF plantations.At the same time, the glpP gene encoding glycerol-3phosphate transporter protein was found to have the highest relative abundance, and it can activate the expression of genes related to the mineralization and transport of glycerol-3-phosphate such as glpF, glpK, glpT, etc. (Darbon et al., 2002).This indicated that the microorganisms in the soil under the three plantations could absorb glycerol-3-phosphate as an alternative phosphorus source to maintain the basic functions of cells, thus coping with phosphorus starvation (Bergkemper et al., 2016).Microbial solubilization of inorganic phosphorus is mainly regulated by the gcd gene encoding quinoprotein glucose dehydrogenase, which is an enzyme that secretes gluconic acid that chelates Ca 2+ , Fe 3+ , and Al 3+ cations, thus dissolving mineral phosphorus that is not directly available to plants and organisms in the soil and releasing available phosphorus (Zeng et al., 2016).The gcd gene not only is ubiquitous in soil but also can be used as a key molecular marker to characterize available phosphorus in soil (Rawat et al., 2020).In this study, gcd in each soil layer under the three stands was dominant, indicating that the microbial solubilization of inorganic phosphorus was the main process providing soil available phosphorus for stands (Liang et al., 2020).At 0-20 cm, gcd was significantly negatively correlated with TP content, which in the PM plantation (0.4 g‧kg −1 ) apparently increased compared with that in the EF plantation (0.28 g‧kg −1 ).Therefore, soil TP is one of the most important factors that significantly increased the relative abundance of gcd after the transformation of the PM plantation into the EF plantation.Previous studies have also pointed out that it is easier for microorganisms to obtain nutrients from the soil with a higher TP content, which leads to a lower investment of microorganisms in the process of phosphorus solubilization, and thus a decrease in relative abundances of related functional genes (Pastore et al., 2020).In addition, the increase in the NO 3 − -N content would increase the demand of microorganisms for inorganic phosphorus, resulting in the improvement of the capacity of microbes to dissolve phosphorus to obtain more inorganic phosphorus for their survival (Xiao et al., 2018).This study found that the relative abundance of the gcd gene was significantly positively correlated with the NO 3 − -N content.The NO 3 − -N content in the 0-20 cm of the PM plantation was significantly lower (0.48 mg‧kg −1 ) than that of the EF plantation (1.56 mg‧kg −1 ), which is another important reason for the higher relative abundance of gcd in the latter.

Effects of stand types on molecular ecological network structure of key genes involved in P transformation
There are many microorganisms associated with complex P-cycling processes in soil (Billah et al., 2019), and molecular ecological networks can reveal the potential interactions of microorganisms (Deng et al., 2012).Compared to simpler networks, in complex networks, microbial communities have more efficient resource utilization and information transmission (Morrien et al., 2017).In this study, the average degree of molecular ecological networks of soil phosphorus cycling microbial functional genes was higher in the PM stand than in the EF stand, but lower than that in the PM/EF stand, indicating that the network complexity was reduced after the transformation of the PM plantation stand into the EF stand but it increased after the transformation into the PM/ EF stand.This is consistent with previous research results, reporting that compared with pure plantations, the molecular ecological network structure of phosphorus cycling functional genes was more complex in the mixed plantation (Cheng et al., 2021).Furthermore, the diversity of tree species, litter, and root exudates in mixed forests was generally higher than in pure forests, which may lead to more complex interactions between soil microorganisms (Haichar et al., 2014).The mean clustering coefficient of molecular networks of soil phosphorus cycling microbial functional genes in the PM stand was lower than that in EF and PM/EF stands, indicating that the level of community organization was the lowest in the PM stand.The lower value of mean path length indicates the closer interaction between genes within this network and the faster diffusion of external interference in networks, resulting in the instability of the network structure (Deng et al., 2012).Higher modularity values indicate a stronger anti-interference ability (Carpenter et al., 2012).The average path length and modularity of the network of the PM stand increased compared to those of the EF stand but they had lower values than those of the PM/EF stand, indicating that soil phosphorus cycling microbial functional genes under this stand type could be more closely functionally related and be more tolerant to environmental disturbances (Yuan et al., 2021), while in EF, the microbial functional genes of soil phosphorus cycling were susceptible to the interference of external environment, and more key genes were needed to maintain the stability of molecular ecological network (Banerjee et al., 2018).In addition, the functional gene networks for the three stands had mainly positive connections, indicating that microbial functional genes involved in soil phosphorus cycling may obtain limited nutrient sources through cooperation (Hoek et al., 2016).
Module hubs and connectors are key genes in functional molecular ecological networks, and they are also critical for maintaining the stability of community structure (Banerjee et al., 2018).In the present study, ugpC was identified as the key gene in the Molecular ecological networks of soil phosphorus cycling microbial functional genes in the PM plantation (A), the EF plantation (B), and the PM/EF plantation (C); Nodes with green, purple, and orange colors represent genes related to P-starvation response regulation, P-solubilization and mineralization, and P-uptake and transport, respectively.The node size indicates the strength of connection with other nodes (genes).Different lines represent the linear correlation between nodes, with solid and dotted lines denoting the positive and negative relations, respectively. 10.3389/fmicb.2024.1419645 Frontiers in Microbiology 14 frontiersin.orgPM plantation, whereas phoU, pstA, ugpE, and gcd were the key genes in the EF plantation, and pqqE and pstB were the key genes in the PM/ EF plantation, all associated with different P cycling processes.
Numerous recent studies demonstrate the potential keystone roles of pyrroloquinoline quinone (pqq) and glucose dehydrogenase (gcd) in phosphorus solubilization, and thus, they are becoming reference biomarkers in research on soil P cycling (Rawat et al., 2020).The pstA, pstB, ugpC, and ugpE genes were expressed by pst and ugp transport systems (pstSCAB and ugpABE, respectively) (Dai et al., 2020).One study has shown that pstB and ugpA are the key functional genes related to soil P cycling (Liu et al., 2023), which conforms to our findings.However, to further verify our results, it should be explored whether the gene coding for the negative regulatory protein phoU plays a key role in P-cycling processes.The modification of key functional genes occurred by the transformation of the PM plantation into EF and PM/EF plantations, which may be due to its influence on the soil microenvironment, thus modulating the response mechanism of microorganisms to environmental disturbances.Overall, this study only collected soil samples in winter, which could not characterize the effects of seasonal changes on the microbial functional genes of soil phosphorus cycling.Soil samples can be collected in different seasons to elucidate the effects of different stand types on soil phosphorus cycling in future.

Conclusion
The composition structure of soil phosphorus cycling microbial functional genes in the PM plantation was significantly different from that in EF and PM/EF plantations in the subtropics of China, mainly due to the difference in soil TP content.After the transformation of the PM plantation stand into the EF stand, the relative abundances of the group of genes involved in inorganic P-solubilization and organic P-mineralization in the topsoil and middle layer significantly increased, and the difference was caused by TP.At the same time, for genes related to P-starvation response regulation, relative abundances in the subsoil increased, and SOC may be the main factor contributing to this difference.The transformation of the PM stand into PM/EF and EF plantation stands had no obvious influence on the dominant phosphorus cycling microbial functional genes in all soil layers of stands.The decrease in the TP content and the increase in the NO 3 − -N content in the topsoil after the transformation of the PM stand into the EF stand are the reasons for the significant increase in the relative abundance of gcd in the EF stand.After the conversion of the PM plantation into the PM/EF plantation, the molecular ecological network structure of soil phosphorus cycling microbial functional genes became more complex, and thus, its stability was enhanced.Moreover, the key genes that play roles in maintaining network stability in different stands exhibited an inconsistent pattern of change.Future studies may be required to focus on the processes of gene transcription and translation, which will help to more comprehensively understand how plantation transformation affects the expression of soil phosphorus cycling microbial functional genes and its potential mechanism.You, Y., Wu, X., Ming, A., Liu, T., Chen, Y., Zhu, H., et al. (2018a).Changes of plant functional group in understory and environmental interpretation in the transformation of typical coniferous plantation to native broadleaved species plantation in south subtropical China.Chin.J. Ecol. 37, 3194-3201. doi: 10.13292/ j.1000-4890.201811.025 You, Y., Xu, H., Wu, X., Zhou, X., Tan, X., Li, M., et al. (2020).Native broadleaf tree species stimulate topsoil nutrient transformation by changing microbial community composition and physiological function, but not biomass in subtropical plantations with low P status.For.Ecol. Manag. 477:118491. doi: 10.1016/j. foreco.2020.118491 Yu, H., Wang, F., Shao, M., Huang, L., Xie, Y., Xu, Y., et al. (2021).Effects of rotations with legume on soil functional microbial communities involved in phosphorus transformation.Front.Microbiol.12:661100.doi: 10.3389/ fmicb.2021.661100Yuan, M. M., Guo, X., Wu, L., Zhang, Y., Xiao, N., Ning, D., et al. (2021).Climate warming enhances microbial network complexity and stability.Nat.Clim. Chang. 11, 343-348. doi: 10.1038/s41558-021-00989-9 Zaidi, A., Khan, M. S., Ahemad, M., andOves, M. (2009).Plant growth promotion by phosphate solubilizing bacteria.Acta Microbiol. Immunol. Hung. 56, 263-284. doi: 10.1556/AMicr.56.2009.3.6 Zeng, Q., Wu, X., and Wen, X. (2016).Effects of soluble phosphate on phosphate-solubilizing characteristics and expression of gcd gene in Pseudomonas frederiksbergensis JW-SD2.Curr.Microbiol. 72,[198][199][200][201][202][203][204][205][206].1007/s00284-015-0938-zZhou, J., Deng, Y., Luo, F., He, Z., Tu, Q., and Zhi, X. (2010).Functional molecular ecological networks.MBio 1, e00169-e00110.doi: 10.1128/mBio.00169-10

FIGURE 2
FIGURE 2Relative abundances of the groups of phosphorus cycling microbial functional genes in different stands at different soil layers; (A) Genes related to P-starvation response regulation; (B) Genes related to P-uptake and transport; and (C) Genes related to P-solubilization and mineralization; PM, Pinus massoniana; EF, Erythrophleum fordii.

FIGURE 4
FIGURE 4The relative abundance of phosphorus cycling microbial functional genes in three studied stands at soil layers of 0-20 (A), 20-40 (B), and 40-60 cm (C); Functional genes with the relative abundance of less than 0.1‰ were eliminated.Different lowercase letters indicate significant differences among different stand types (p < 0.05).

FIGURE 7
FIGURE 7Topological clustering of soil phosphorus cycling microbial functional genes in the PM plantation (A), the EF plantation (B), and the PM/EF plantation (C).

TABLE 1
The general information of the three studied planted forests.

TABLE 2
Classification of soil phosphorus cycling microbial functional genes.

TABLE 3
Comparison of soil physicochemical properties among three different plantations (mean ± SD, n = 3).-N, Ammonium nitrogen; and AP, Available phosphorus.Different lowercase letters indicate significant differences among different stand types for the same soil layer, whereas different uppercase letters denote significant differences among different soil layers for the same stand type at the 0.05 significance level.

TABLE 4
Results of the two-way ANOVA showing the impacts of stand types and soil layers on phosphorus cycling microbial functional genes.

TABLE 5
Topological properties of molecular ecological networks of phosphorus cycling microbial functional genes in three studied stands.

TABLE 6
Module hubs and connecters of the soil phosphorus cycling microbial functional genes molecular ecological networks in different stands.